$title balance_benchmark
* filters small values from the benchmark dataset and rebalances the sam



* ------------------------------------------------------------------------------
* set options
* ------------------------------------------------------------------------------

* set default benchmark file
$if not set benchmark_file $set benchmark_file "data\default_aggregation.gdx"

* set default output file
$if not set output_file $set output_file "data\default_aggregation.gdx"

* set default for disaggregating activities based on extant capital
$if not set filter_small $set filter_small 1

* set default for taxes to change in the balancing
$if not set include_taxes $set include_taxes 0

* smallest value allowed in benchmark
$if not set threshold $set threshold 5e-4

* set default to minimize fractional deviations not absolute deviations
$if not set frac_deviations $set frac_deviations 1

* set the default option to use the balanced growth investment constraint
$if not set balanced_growth $set balanced_growth 1

* set default growth rate for the balanced growth constraint
$if not set growth_rate $set growth_rate 0.02

* set default depreciation rate for the balanced growth constraint
$if not set delta $set delta 0.05

* set default real interest rate for the balanced growth constraint
$if not set rbar $set rbar 0.05

* set default interval for the model timestep in years
$if not set interval $set interval 5



* ------------------------------------------------------------------------------
* load model sets
* ------------------------------------------------------------------------------

* basic model sets
sets
  s                     sectors
  r                     regions
  h                     households
  trd                   trade markets
  f                     primary factors
  t                     time periods;

* set aliases
alias
  (s,ss,sss), (h,hh);

$gdxin %benchmark_file%
$load r s h trd f t



* ------------------------------------------------------------------------------
* load parameters
* ------------------------------------------------------------------------------

* benchmark data
parameters
    y0(r,s)         benchmark - output
    ld0(r,s)        benchmark - labor demand
    kd0(r,s)        benchmark - capital demand
    id0(r,s,ss)     benchmark - intermediate demand
    d0(r,s)         benchmark - domestic supply
    m0(r,s,trd)     benchmark - imports
    x0(r,s,trd)     benchmark - exports
    c0(r,h)         benchmark - consumption bundle
    cd0(r,s,h)      benchmark - consumer demand
    i0(r,s)         benchmark - investment demand
    invh0(r,h)      benchmark - household investment
    n0(r,h)         benchmark - households
    g0(r,s)         benchmark - government demand
    le0(r,h)        benchmark - household labor supply
    re0(r,h)        benchmark - capital rental income
    bopdef(r,h)     benchmark - international balance of payments
    incadj0(r,h)    benchmark - government to household lump sum transfer
    ty0(r,s)        tax rate - production
    tl0(r)          tax rate - labor use
    tk0(r)          tax rate - capital use
    tc0(r)          tax rate - consumption;

* Load benchmark data
$load id0 ld0 kd0 cd0 x0 m0 i0 g0 re0 le0 invh0 n0 ty0 tl0 tc0 tk0

* calculate remaining variables from the original benchmark
y0(r,s)      = (sum(ss, id0(r,ss,s))+(1+tl0(r))*ld0(r,s)+(1+tk0(r))*kd0(r,s))/(1-ty0(r,s));
c0(r,h)      = sum(s, cd0(r,s,h))*(1+tc0(r));
bopdef(r,h)  = (sum(s, m0(r,s,"ftrd")-x0(r,s,"ftrd")))*c0(r,h)/sum(hh, c0(r,hh));
incadj0(r,h) = c0(r,h)+invh0(r,h)-re0(r,h)-le0(r,h)-bopdef(r,h);

* skip to the balancing routine if filtering small values not requested
$if not %filter_small%==1 $goto balance



* ------------------------------------------------------------------------------
* filter out small values
* ------------------------------------------------------------------------------

parameters
  non_zeros_before      non-zero benchmark elements before filter
  non_zeros_after       non-zero benchmark elements after filter
  reduction             fraction of non-zeros elements filtered;

* non-zero benchmark elements before filter
non_zeros_before = sum((r,s,h,trd,ss), 1$ld0(r,s)
                                      +1$kd0(r,s)
                                      +1$id0(r,s,ss)
                                      +1$cd0(r,s,h)
                                      +1$x0(r,s,trd)
                                      +1$m0(r,s,trd)
                                      +1$i0(r,s)
                                      +1$g0(r,s)
                                      +1$re0(r,h)
                                      +1$le0(r,h)
                                      +1$invh0(r,h)
                                      +1$ty0(r,s)
                                      +1$tl0(r)
                                      +1$tc0(r)
                                      +1$tk0(r));

* filter out small output values
y0(r,s)$(y0(r,s)<%threshold%) = 0.0;

* there is no intermediate or factor demand or exports if there is no output
id0(r,ss,s)$(y0(r,s) eq 0) = 0.0;
ld0(r,s)$(y0(r,s) eq 0)    = 0.0;
kd0(r,s)$(y0(r,s) eq 0)    = 0.0;
x0(r,s,trd)$(y0(r,s) eq 0) = 0.0;

* filter out other small values
i0(r,s)$(i0(r,s)<%threshold%) = 0.0;
g0(r,s)$(g0(r,s)<%threshold%) = 0.0;
cd0(r,s,h)$(cd0(r,s,h)<%threshold%) = 0.0;
m0(r,s,trd)$(m0(r,s,trd)<%threshold%) = 0.0;
x0(r,s,trd)$(x0(r,s,trd)<%threshold%) = 0.0;

* filter out inputs that are a small share of costs
id0(r,ss,s)$((id0(r,ss,s)/y0(r,s))<%threshold%) = 0.0;

* non-zero benchmark elements after filter
non_zeros_after = sum((r,s,h,trd,ss), 1$ld0(r,s)
                                     +1$kd0(r,s)
                                     +1$id0(r,s,ss)
                                     +1$cd0(r,s,h)
                                     +1$x0(r,s,trd)
                                     +1$m0(r,s,trd)
                                     +1$i0(r,s)
                                     +1$g0(r,s)
                                     +1$re0(r,h)
                                     +1$le0(r,h)
                                     +1$invh0(r,h)
                                     +1$ty0(r,s)
                                     +1$tl0(r)
                                     +1$tc0(r)
                                     +1$tk0(r));

*fraction of non-zeros elements filtered
reduction = non_zeros_after/non_zeros_before-1;

display non_zeros_before, non_zeros_after, reduction;

$label balance



* ------------------------------------------------------------------------------
* balance benchmark dataset
* ------------------------------------------------------------------------------

* rebalanced benchmark data
positive variables
  y1(r,s)               benchmark - output
  id1(r,s,ss)           benchmark - intermediate demand
  ld1(r,s)              benchmark - labor demand
  kd1(r,s)              benchmark - capital demand
  m1(r,s,trd)           benchmark - imports
  x1(r,s,trd)           benchmark - exports
  cd1(r,s,h)            benchmark - consumer demand
  i1(r,s)               benchmark - investment demand
  invh1(r,h)            benchmark - household investment
  g1(r,s)               benchmark - government demand
  le1(r,h)              benchmark - household labor supply
  re1(r,h)              benchmark - capital rental income
  tl1(r)                tax rate - labor used
  tk1(r)                tax rate - capital use
  tc1(r)                tax rate - consumption
  ty1(r,s)              tax rate - production;

* endogenous economic account variables
variables
  incadj1(r,h)          benchmark - government to household lump sum transfer;

* objective function
variables
  ssd                   sum of the squared deviation from starting becnhmark;

* objective function and constraints
equations
  sseq                  sum of the squared deviation from starting becnhmark
  ldeq                  labor market clearnace condition
  kdeq                  capital market clearnace condition
  ideq                  investment market clearnace condition
  mceq                  market clearance conditions
  zpeq                  zero profit conditions
  dteq                  inter-regional trade balance constraints
  bceq                  household budget constraints
  gbeq                  government budget constraint
  dmeq                  weakly positive domestic own use
  bdeq                  maintain balance of payments household share out
  inv_ss                steady state investment level constraint [optional];

* set starting values to the original benchmark
y1.l(r,s)      = y0(r,s);
id1.l(r,s,ss)  = id0(r,s,ss);
ld1.l(r,s)     = ld0(r,s);
kd1.l(r,s)     = kd0(r,s);
cd1.l(r,s,h)   = cd0(r,s,h);
x1.l(r,s,trd)  = x0(r,s,trd);
m1.l(r,s,trd)  = m0(r,s,trd);
i1.l(r,s)      = i0(r,s);
g1.l(r,s)      = g0(r,s);
le1.l(r,h)     = le0(r,h);
re1.l(r,h)     = re0(r,h);
invh1.l(r,h)   = invh0(r,h);
incadj1.l(r,h) = incadj0(r,h);
tl1.l(r)       = tl0(r);
tk1.l(r)       = tk0(r);
tc1.l(r)       = tc0(r);
ty1.l(r,s)     = ty0(r,s);

* maintain original sparsity
y1.fx(r,s)$(y0(r,s) eq 0)            = 0;
id1.fx(r,s,ss)$(id0(r,s,ss) eq 0)    = 0;
ld1.fx(r,s)$(ld0(r,s) eq 0)          = 0;
kd1.fx(r,s)$(kd0(r,s) eq 0)          = 0;
cd1.fx(r,s,h)$(cd0(r,s,h) eq 0)      = 0;
x1.fx(r,s,trd)$(x0(r,s,trd) eq 0)    = 0;
m1.fx(r,s,trd)$(m0(r,s,trd) eq 0)    = 0;
i1.fx(r,s)$(i0(r,s) eq 0)            = 0;
g1.fx(r,s)$(g0(r,s) eq 0)            = 0;
le1.fx(r,h)$(le0(r,h) eq 0)          = 0;
re1.fx(r,h)$(re0(r,h) eq 0)          = 0;
invh1.fx(r,h)$(invh0(r,h) eq 0)      = 0;
incadj1.fx(r,h)$(incadj0(r,h) eq 0)  = 0;
tl1.fx(r)$(tl0(r) eq 0)              = 0;
tk1.fx(r)$(tk0(r) eq 0)              = 0;
tc1.fx(r)$(tc0(r) eq 0)              = 0;
ty1.fx(r,s)$(ty0(r,s) eq 0)          = 0;

* fix taxes if requested
if(%include_taxes% ne 1,
  tl1.fx(r)   = tl0(r);
  tk1.fx(r)   = tk0(r);
  tc1.fx(r)   = tc0(r);
  ty1.fx(r,s) = ty0(r,s);
);

$if not %frac_deviations%==1 $goto absolute_deviations

* sum of the squared percentage deviation from starting becnhmark
sseq..  ssd =e=  sum((r,s,ss)$id0(r,s,ss),  abs(id0(r,s,ss)) *power(id1(r,s,ss)  /id0(r,s,ss) -1,2))
                +sum((r,s)$ld0(r,s),        abs(ld0(r,s))    *power(ld1(r,s)     /ld0(r,s)    -1,2))
                +sum((r,s)$kd0(r,s),        abs(kd0(r,s))    *power(kd1(r,s)     /kd0(r,s)    -1,2))
                +sum((r,s,trd)$m0(r,s,trd), abs(m0(r,s,trd)) *power(m1(r,s,trd)  /m0(r,s,trd) -1,2))
                +sum((r,s,trd)$x0(r,s,trd), abs(x0(r,s,trd)) *power(x1(r,s,trd)  /x0(r,s,trd) -1,2))
                +sum((r,s,h)$cd0(r,s,h),    abs(cd0(r,s,h))  *power(cd1(r,s,h)   /cd0(r,s,h)  -1,2))
                +sum((r,s)$g0(r,s),         abs(g0(r,s))     *power(g1(r,s)      /g0(r,s)     -1,2))
                +sum((r,s)$i0(r,s),         abs(i0(r,s))     *power(i1(r,s)      /i0(r,s)     -1,2))
                +sum((r,h)$le0(r,h),        abs(le0(r,h))    *power(le1(r,h)     /le0(r,h)    -1,2))
                +sum((r,h)$re0(r,h),        abs(re0(r,h))    *power(re1(r,h)     /re0(r,h)    -1,2))
                +sum((r,h)$invh0(r,h),      abs(invh0(r,h))  *power(invh1(r,h)   /invh0(r,h)  -1,2))
                +sum((r,h)$incadj0(r,h),    abs(incadj0(r,h))*power(incadj1(r,h) /incadj0(r,h)-1,2))
                +sum((r)$tl0(r),                              power(tl1(r)       /tl0(r)      -1,2))
                +sum((r)$tk0(r),                              power(tk1(r)       /tk0(r)      -1,2))
                +sum((r)$tc0(r),                              power(tc1(r)       /tc0(r)      -1,2))
                +sum((r,s)$ty0(r,s),                          power(ty1(r,s)     /ty0(r,s)    -1,2));

$goto constraints
$label absolute_deviations

* sum of the squared deviation from starting becnhmark
sseq..  ssd =e=  sum((r,s,ss),  power(id0(r,s,ss)  -id1(r,s,ss) ,2))
                +sum((r,s),     power(ld0(r,s)     -ld1(r,s)    ,2))
                +sum((r,s),     power(kd0(r,s)     -kd1(r,s)    ,2))
                +sum((r,s,trd), power(m0(r,s,trd)  -m1(r,s,trd) ,2))
                +sum((r,s,trd), power(x0(r,s,trd)  -x1(r,s,trd) ,2))
                +sum((r,s,h),   power(cd0(r,s,h)   -cd1(r,s,h)  ,2))
                +sum((r,s),     power(g0(r,s)      -g1(r,s)     ,2))
                +sum((r,s),     power(i0(r,s)      -i1(r,s)     ,2))
                +sum((r,h),     power(le0(r,h)     -le1(r,h)    ,2))
                +sum((r,h),     power(re0(r,h)     -re1(r,h)    ,2))
                +sum((r,h),     power(invh0(r,h)   -invh1(r,h)  ,2))
                +sum((r,h),     power(incadj0(r,h) -incadj1(r,h),2));

$label constraints

* regional factor endowments equal factor demands
ldeq(r)..
  sum(h, le1(r,h)) =e= sum(s, ld1(r,s));
kdeq(r)..
  sum(h, re1(r,h)) =e= sum(s, kd1(r,s));

* ensure that regional household investment equals investment demands
ideq(r)..
  sum(h, invh1(r,h)) =e= sum(s, i1(r,s));

* market clearance conditions
mceq(r,s)..
  y1(r,s)+sum(trd, m1(r,s,trd)-x1(r,s,trd)) =e= sum(ss, id1(r,s,ss))+g1(r,s)
                                                +i1(r,s)+sum(h, cd1(r,s,h));

* zero profit conditions
zpeq(r,s)..
  (1-ty1(r,s))*y1(r,s) =e= sum(ss, id1(r,ss,s))+ld1(r,s)*(1+tl1(r))
                           +kd1(r,s)*(1+tk1(r));

* inter-regional trade balance constraints
dteq(s)..
  sum(r, m1(r,s,"dtrd")) =e= sum(r, x1(r,s,"dtrd"));

* household budget constraints
bceq(r,h)..
  sum(s, (1+tc1(r))*cd1(r,s,h))+invh1(r,h) =e= re1(r,h)+le1(r,h)+incadj1(r,h)
                                               +bopdef(r,h);

* government budget constraint
gbeq..
  sum((r,s), g1(r,s))+sum((r,h), incadj1(r,h))
    =e= sum((r,s), tc1(r)*sum(h, cd1(r,s,h))+tk1(r)*kd1(r,s)+tl1(r)*ld1(r,s)
        +ty1(r,s)*y1(r,s));

* (weakly) positive domestic own use
dmeq(r,s)..
  0.99*y1(r,s) =g= sum(trd, x1(r,s,trd));

* maintain balance of payments household share out
bdeq(r,h)..
  bopdef(r,h) =e= (sum(s, m1(r,s,"ftrd")-x1(r,s,"ftrd")))
                  *sum(s,cd1(r,s,h))/sum((s,hh), cd1(r,s,hh));

* if a balanced growth path is requested constrain investment to be consistent
* with that assumption, if not instert a placebo condition for the constraint
* this seems easier then adding two differnt model definitions
$if not %balanced_growth%==1 $goto generic_inv_ss

* steady-state investment level
inv_ss(r)..
  sum(s, i1(r,s)) =e= ((1+%growth_rate%)**%interval%-1+1-(1-%delta%)**%interval%)
                      *sum(h, re1(r,h))/((1+%rbar%)**%interval%-1+1-(1-%delta%)**%interval%);

$goto done_inv_ss
$label generic_inv_ss

* placebo constraint
inv_ss(r)..
  1 =e= 1;

$label done_inv_ss

* rebalance the data
option nlp=conopt;
model rebalance SAM rebalance / ALL /;
solve rebalance using nlp minimizing ssd;

* intermediates and results for rebalance statistics
parameter
  rebalStats(*,*)       statistics regarding the rebalance changes
  idr(r,s,ss)           benchmark - intermediate demand
  ldr(r,s)              benchmark - labor demand
  kdr(r,s)              benchmark - capital demand
  mr(r,s,trd)           benchmark - imports
  xr(r,s,trd)           benchmark - exports
  cdr(r,s,h)            benchmark - consumer demand
  ir(r,s)               benchmark - investment demand
  invhr(r,h)            benchmark - household investment
  gr(r,s)               benchmark - government demand
  ler(r,h)              benchmark - household labor supply
  rer(r,h)              benchmark - capital rental income
  incadjr(r,h)          benchmark - government to household lump sum transfer
  tlr(r)                tax rate - labor
  tkr(r)                tax rate - capital
  tcr(r)                tax rate - consumption
  tyr(r,s)              tax rate - production;

* default the ratios to one - no change
idr(r,s,ss)  = 1;
ldr(r,s)     = 1;
kdr(r,s)     = 1;
cdr(r,s,h)   = 1;
xr(r,s,trd)  = 1;
mr(r,s,trd)  = 1;
ir(r,s)      = 1;
gr(r,s)      = 1;
ler(r,h)     = 1;
rer(r,h)     = 1;
invhr(r,h)   = 1;
incadjr(r,h) = 1;
tlr(r)       = 1;
tkr(r)       = 1;
tcr(r)       = 1;
tyr(r,s)     = 1;

* ratio of the individual parameter changes
idr(r,s,ss)$id0(r,s,ss)   = id1.l(r,s,ss)  / id0(r,s,ss);
ldr(r,s)$ld0(r,s)         = ld1.l(r,s)     / ld0(r,s);
kdr(r,s)$kd0(r,s)         = kd1.l(r,s)     / kd0(r,s);
cdr(r,s,h)$cd0(r,s,h)     = cd1.l(r,s,h)   / cd0(r,s,h);
xr(r,s,trd)$x0(r,s,trd)   = x1.l(r,s,trd)  / x0(r,s,trd);
mr(r,s,trd)$m0(r,s,trd)   = m1.l(r,s,trd)  / m0(r,s,trd);
ir(r,s)$i0(r,s)           = i1.l(r,s)      / i0(r,s);
gr(r,s)$g0(r,s)           = g1.l(r,s)      / g0(r,s);
ler(r,h)$le0(r,h)         = le1.l(r,h)     / le0(r,h);
rer(r,h)$re0(r,h)         = re1.l(r,h)     / re0(r,h);
invhr(r,h)$invh0(r,h)     = invh1.l(r,h)   / invh0(r,h);
incadjr(r,h)$incadj0(r,h) = incadj1.l(r,h) / incadj0(r,h);
tlr(r)$tl0(r)             = tl1.l(r)       / tl0(r);
tkr(r)$tk0(r)             = tk1.l(r)       / tk0(r);
tcr(r)$tc0(r)             = tc1.l(r)       / tc0(r);
tyr(r,s)$ty0(r,s)         = ty1.l(r,s)     / ty0(r,s);

* maximum pre-post rebalance ratios
rebalStats("id","ratio_max")     = smax((r,s,ss),  idr(r,s,ss)  );
rebalStats("ld","ratio_max")     = smax((r,s),     ldr(r,s)     );
rebalStats("kd","ratio_max")     = smax((r,s),     kdr(r,s)     );
rebalStats("cd","ratio_max")     = smax((r,s,h),   cdr(r,s,h)   );
rebalStats("x","ratio_max")      = smax((r,s,trd), xr(r,s,trd)  );
rebalStats("m","ratio_max")      = smax((r,s,trd), mr(r,s,trd)  );
rebalStats("i","ratio_max")      = smax((r,s),     ir(r,s)      );
rebalStats("g","ratio_max")      = smax((r,s),     gr(r,s)      );
rebalStats("le","ratio_max")     = smax((r,h),     ler(r,h)     );
rebalStats("ke","ratio_max")     = smax((r,h),     rer(r,h)     );
rebalStats("invh","ratio_max")   = smax((r,h),     invhr(r,h)   );
rebalStats("incadj","ratio_max") = smax((r,h),     incadjr(r,h) );
rebalstats("tl","ratio_max")     = smax((r),       tlr(r)       );
rebalstats("tk","ratio_max")     = smax((r),       tkr(r)       );
rebalstats("tc","ratio_max")     = smax((r),       tcr(r)       );
rebalstats("ty","ratio_max")     = smax((r,s),     tyr(r,s)     );

* minimum pre-post rebalance ratios
rebalStats("id","ratio_min")     = smin((r,s,ss),  idr(r,s,ss)  );
rebalStats("ld","ratio_min")     = smin((r,s),     ldr(r,s)     );
rebalStats("kd","ratio_min")     = smin((r,s),     kdr(r,s)     );
rebalStats("cd","ratio_min")     = smin((r,s,h),   cdr(r,s,h)   );
rebalStats("x","ratio_min")      = smin((r,s,trd), xr(r,s,trd)  );
rebalStats("m","ratio_min")      = smin((r,s,trd), mr(r,s,trd)  );
rebalStats("i","ratio_min")      = smin((r,s),     ir(r,s)      );
rebalStats("g","ratio_min")      = smin((r,s),     gr(r,s)      );
rebalStats("le","ratio_min")     = smin((r,h),     ler(r,h)     );
rebalStats("ke","ratio_min")     = smin((r,h),     rer(r,h)     );
rebalStats("invh","ratio_min")   = smin((r,h),     invhr(r,h)   );
rebalStats("incadj","ratio_min") = smin((r,h),     incadjr(r,h) );
rebalstats("tl","ratio_min")     = smin((r),       tlr(r)       );
rebalstats("tk","ratio_min")     = smin((r),       tkr(r)       );
rebalstats("tc","ratio_min")     = smin((r),       tcr(r)       );
rebalstats("ty","ratio_min")     = smin((r,s),     tyr(r,s)     );

* absolute deviations
idr(r,s,ss)  = abs(id1.l(r,s,ss)  - id0(r,s,ss)  );
ldr(r,s)     = abs(ld1.l(r,s)     - ld0(r,s)     );
kdr(r,s)     = abs(kd1.l(r,s)     - kd0(r,s)     );
cdr(r,s,h)   = abs(cd1.l(r,s,h)   - cd0(r,s,h)   );
xr(r,s,trd)  = abs(x1.l(r,s,trd)  - x0(r,s,trd)  );
mr(r,s,trd)  = abs(m1.l(r,s,trd)  - m0(r,s,trd)  );
ir(r,s)      = abs(i1.l(r,s)      - i0(r,s)      );
gr(r,s)      = abs(g1.l(r,s)      - g0(r,s)      );
ler(r,h)     = abs(le1.l(r,h)     - le0(r,h)     );
rer(r,h)     = abs(re1.l(r,h)     - re0(r,h)     );
invhr(r,h)   = abs(invh1.l(r,h)   - invh0(r,h)   );
incadjr(r,h) = abs(incadj1.l(r,h) - incadj0(r,h) );
tlr(r)       = abs(tl1.l(r)       - tl0(r)       );
tkr(r)       = abs(tk1.l(r)       - tk0(r)       );
tcr(r)       = abs(tc1.l(r)       - tc0(r)       );
tyr(r,s)     = abs(ty1.l(r,s)     - ty0(r,s)     );

* max absolute deviations
rebalStats("id","dev_max")     = smax((r,s,ss),  idr(r,s,ss)  );
rebalStats("ld","dev_max")     = smax((r,s),     ldr(r,s)     );
rebalStats("kd","dev_max")     = smax((r,s),     kdr(r,s)     );
rebalStats("cd","dev_max")     = smax((r,s,h),   cdr(r,s,h)   );
rebalStats("x","dev_max")      = smax((r,s,trd), xr(r,s,trd)  );
rebalStats("m","dev_max")      = smax((r,s,trd), mr(r,s,trd)  );
rebalStats("i","dev_max")      = smax((r,s),     ir(r,s)      );
rebalStats("g","dev_max")      = smax((r,s),     gr(r,s)      );
rebalStats("le","dev_max")     = smax((r,h),     ler(r,h)     );
rebalStats("ke","dev_max")     = smax((r,h),     rer(r,h)     );
rebalStats("invh","dev_max")   = smax((r,h),     invhr(r,h)   );
rebalStats("incadj","dev_max") = smax((r,h),     incadjr(r,h) );
rebalstats("tl","dev_max")     = smax((r),       tlr(r)       );
rebalstats("tk","dev_max")     = smax((r),       tkr(r)       );
rebalstats("tc","dev_max")     = smax((r),       tcr(r)       );
rebalstats("ty","dev_max")     = smax((r,s),     tyr(r,s)     );

* display the statistics for the rebalance
display rebalStats;

* save the updated parameters
id0(r,s,ss) = id1.l(r,s,ss);
ld0(r,s)    = ld1.l(r,s);
kd0(r,s)    = kd1.l(r,s);
cd0(r,s,h)  = cd1.l(r,s,h);
x0(r,s,trd) = x1.l(r,s,trd);
m0(r,s,trd) = m1.l(r,s,trd);
i0(r,s)     = i1.l(r,s);
g0(r,s)     = g1.l(r,s);
le0(r,h)    = le1.l(r,h);
re0(r,h)    = re1.l(r,h);
invh0(r,h)  = invh1.l(r,h);
tl0(r)      = tl1.l(r);
tk0(r)      = tk1.l(r);
tc0(r)      = tc1.l(r);
ty0(r,s)    = ty1.l(r,s);

* if a balanced growth path based benchmark was requested
* regional investment is shared out to households by assuming its consistent
* with a balanced growth path with government transfers balancing the budget
* constraint following goulder and hafstead (2018) appendix d
$if not %balanced_growth%==1 $goto write_benchmark
invh0(r,h) = ((1+%growth_rate%)**%interval%-1+1-(1-%delta%)**%interval%)
             *re0(r,h)/((1+%rbar%)**%interval%-1+1-(1-%delta%)**%interval%);
$label write_benchmark



* ------------------------------------------------------------------------------
* save benchmark dataset
* ------------------------------------------------------------------------------

execute_unload "%output_file%", f,s,r,h,t,id0,ld0,kd0,cd0,x0,
                                m0,i0,g0,re0,le0,invh0,n0
                                ty0,tl0,tc0,tk0,trd;
